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C/3 , Abstract. The electronic structure calculations based upon energy density functionals are 

highly successful and widely used both in solid state physics and quantum chemistry. Moreover, 
the Hohenberg-Kohn theorems and the Kohn-Sham method provide them with a firm basis. 
However, several basic issues are not solved, and hamper the progress to achieve high accuracy. 
In this paper we focus on the conceptual problem of size consistency, basing our analysis on the 
^ , non-intensive character of the (spin) electronic density in the presence of degeneracy. We also 

briefly discuss some of the issues concerning fractional electron numbers from the same point 
of view, analyzing the behavior of the exact functionals for the He and the Hooke's atom series 
, when the number of electrons fluctuates between one and two. 
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1. Introduction 

^ I Density functional tlieory (see, e.g., pQ) (DFT) is by now the most popular method for electronic 

structure calculations in condensed matter physics and quantum chemistry, because of its unique 
combination of low computational cost and reasonable accuracy for many molecules and solids. 
However, despite its large success in scientific areas ranging from material science to biology, 
psj . several basic issues in DFT are still unsolved, and hamper futher developments towards high 

C3 I accuracy. Besides, some of these issues are often ovelooked or simply ignored. We believe that 

raising them is a necessary step towards further improvement of DFT performances. To this 
purpose, in this paper we concentrate on a critical issue in this regard, namely the problem of 
size consistency in DFT in the presence of degeneracy. 
^ ■ The Hohenberg-Kohn theorems [2] and the Kohn-Sham method [3] provide a firm basis for 

^ I DFT calculations, in which the ground state energy and density of a many-electron system is 

obtained by using energy density functionals that can be rigorously defined. Unfortunately, the 
exact definition does not give a prescription that can be followed in practice. Approximations 
are needed for the exchange-correlation energy -E'xci'^] as a functional of the electronic density 
n(r), together with the corresponding exchange-correlation part of the Kohn-Sham potential, 
i.e., the functional derivative of -E'xc[?T'] with respect to n(r). It is worth stressing that none of 
the available approximate i?xc[?T-] satisfies two basic requirements: 

• to provide reasonable error estimates; 

• to allow systematic improvement, meaning that one knows how to reduce the errors by 
means of a well defined procedure. 



Moreover, size consistency in DFT, which is often taken for granted (at least for systems 
composed of closed-shell fragments) can still be an issue, as discussed in the following sections. 
This paper is organized as follows. Section [2] discusses the general form of many approximate 
energy density functionals, focusing on the size consistency problem in the Hohenberg-Kohn 
framework. A similar analysis within the Kohn-Sham theory is then carried out in Sec. [3l 
Systems with fractional electron number are discussed from the same point of view in Sec. [H 
The last Sec. \E\ is devoted to concluding remarks. 



2. Size consistency and density functional theory 

The Hohenberg-Kohn theorems [2] state that the energy of a many-electron system is a 
variational functional of the density, E[n], given by the sum of a universal part, F[n], and a 
linear functional determined by the external potential Vne = ^iVnei^^i) acting on the electrons, 

E[n] = F[n] + j dr t>„e(r) n(r). (1) 

Although the great success of DFT comes essentially from the introduction of the Kohn-Sham 
kinetic energy functional, orbital-free DFT is appealing for its lower computational cost, and is 
nowadays an active field of research (see, e.g., [Ij). It is thus worth to first analyze the size- 
consistency issue in the pure Hohenberg-Kohn framework. In orbital-free DFT the universal 
functional F[n] of Eq. ([1]) is rewritten as 

F[n] = E^^,[n] + Eu[n], (2) 

where i?H[f^] is the usual Hartree classical repulsion energy, Eii[n] = ^ f dr f dr' n{r)n{r')\r — 
r'l^^, and E'kxci?^] is the remaining part of the energy, the kinetic and exchange correlation 
functional, that needs to be approximated. The simplest approximations to -E'kxci?^'] have the 
form 

E^If''[n]= I dvf{n(v),\Vn{v)\,...), (3) 



where the function / is chosen to yield accurate properties for selected systems, e.g., the energy 
of the uniform electron gas, and/or to satisfy some known exact constraints, like scaling relations. 

The requirement of size consistency for a quantum chemistry method is that the result for 
the energy E{A + B) of two non-interacting systems A and B (e.g., which are at infinite distance 
from each other) be equal to the sum of their individual energies, E{A) + E{B). Size consistency 
is evidently crucial when computing dissociation energies. 

When approximations of the form of Eq. ^ are used, one is tempted to believe that size- 
consistency is guaranteed if the function / is intensive, i.e., if its value in the domain of space 
pertaining to system A, 0^, is not changed by the presence of the system B very far from A. A 
corresponding statement can be made for the integration VLb over the region of system B. As 
the integral in the composite system is the sum over the regions of the individual (sub-)systems, 

/dr/(n(r),|Vn(r)l,...) = / dvf{n{v),\Vn{v)\,...)+f dr/(n(r), lVn(r)|, ...) (4) 

size consistency would be guaranteed (for the Hartree and the external potential terms the 
same partition obviously holds for non-interacting subsystems). The basic belief behind this 
statement is that the density itself be an intensive quantity, i.e., that the total density nA+B{T^) 
of the composite system A + B he equal to 

, . / nA(r) r e 



However, this is not true when one of the two systems (or both) has a degenerate ground state 
with different densities. In this case, even an infinitesimal interaction with the other system can 
select one of the states (or a specific ensemble of some of the degenerate states). This change is 
not infinitesimal, and depends on the nature of the other system. Simple examples are diatomic 
molecules in which the individual atoms have partially filled shells with i > (e.g. B2, C2, 
...). An even simpler example is one of such atoms (with degenerate non-spherical densities) 
perturbed by a proton placed at a large distance. Thus, for degenerate systems, Eq. 1^ should 
be replaced by 

where nAi{r) and nBi{'r) are, respectively, the i"^ degenerate densities of A and B, and the 
notation Wi{B) (and Wi{A)) explicits the dependence of the weights of each ensemble (^- Wi = \ 
and < ifj < 1) on the presence of the other system, even if very far away. Equation ([6]) shows 
the non intensive character of the density in the presence of degeneracy. 

The exact F[n]^ of course, would preserve the degeneracy so that it would give for any linear 
combination of the degenerate densities, say, nAii^) the same energy for the system A, leading 
to size consistency. But none of the available approximate functionals is able to preserve the 
degeneracy of the physical system, as none is invariant within the set of degenerate densities. 
Typically, when treating the isolated systems with an approximate functional, one obtains a 
lower energy with one of the degenerate densities, or for a particular ensemble. This means 
that the approximation is size consistent only for some very specific choices of A and B, not in 
general. 

Notice that here we consider the simple case in which different degenerate wavefunctions yield 
different densities. The interesting case of different denegerate wavefunctions corresponding to 
the same density is deeply analyzed in [5]. 

3. Kohn-Sham framework 

In their foundational work, Kohn and Sham j3j split the functional -Ekxc M of Eq. ^ into 

SkxcN =r,[n] + s,cN, (7) 

where Ts[r}\ is the kinetic energy of a system of non-interacting fermions with density n(r) 

T^fn] = max |nnn($|r-hV"|'l>) -ydrn(r)t;(r)| , with F = ^ t;(rj), (8) 

where ^ is in most cases a single Slater determinant. The spin-orbitals (j)i entering in $ are 
determined via the self-consistent equations 



-^V^ -h fH(r)[n] + ■Uxc(r)[n] + Une(r) 



0i(r) = ei(/)i(r), n(r) = ^ /j|(/)i(r)p, 



The Kohn-Sham potential wks = V}i[n\ + v^dn] + Vne is the maximizing potential of Eq. dHj). 
Common approximations for the exchange-correlation energy E^c [t^] are typically of the form of 
Eq. ([3]), so that as far as size consistency in the presence of degeneracy is concerned, we can 
apply to the Kohn-Sham -E'xci'^] the same considerations of the previous section, although in 
this case the situation is complicated by the presence of the functional T^fra], which depends 



on the density in a rather complex way. Ts[n] can have different values for different densities 
ni{r) that are degenerate in the physical system. The work of Fertig and Kohn [7] clearly shows 
that to different degenerate nj(r) can correspond different Kohn-Sham potentials UKS,i(r). For 
example, in an open shell atom with several degenerate non-spherical densities we have different 
non-spherical Kohn-Sham potentials, one for each symmetry. And for an ensemble of these 
densities, we usually need yet another Kohn-Sham potential [8]. 

The requirement for approximate E^c [n] to recover the degeneracy of the physical system has 
been recognized by several authors. Usually, it is written in the form (see, e.g., Ref. [S]) 

Tsl^iWi Hi] + E^ci'^iWi Ui] + E-nlT^iWi Ui] = T^iWi (T^n^] + -Exc [""-«] + Eu[ni]) , (10) 

where rii are the densities of the physical system corresponding to a set of orthonormalized 
degenerate ground-state wavefunctions , and Eq. (fTOl) should hold for any set of the weights 
Wi. However, as noted in Refs. (IHIIS], imposing Eq. (fTO]l to approximate Exc[n] is probably a 
daunting task: the Hartree energy E'nfSjWj rij] contains cross terms ij, and must be compensated 
by a complex interplay between Tg and Ey^c- This is illustrated with simple examples in the next 
SecH 

In practical Kohn-Sham calculations, the individual spin densities, n|(r) and ^^J,(^), are used 
as two indipendent variables. In this case the total energy is rewritten as 

E[n^,ni]=Ts[n^,ni]+E-n[n]+E^c[n^,ni]+ J drv{r)n{r). (11) 

As the total density, the spin densities are also non-intensive. Besides, degeneracy can occur 
more often than when considering the total density only. The simplest example is the hydrogen 
atom, which has two degenerate set of spin densities, = n, = and = 0, = n. In 
the streched hydrogen molecule, we have the equi-ensemble of the two on each atom. 

4. Fractional number of electrons 

The behavior of density functionals in the presence of a fractional number of electrons has been 
widely investigated [HI 121 ES] • Nowadays, there is a renovated interest in this issue (see, e.g., 
|14t [TCI [TBj ) , which has lead to the definition of the "many-electron self- interaction error" . 

Fractional number of electrons can be recast in the same class of problems (degeneracy and 
size-consistency) discussed in the previous sections (see also Ref. [9J). To understand why, 
consider the simple example of the stretched Hcg" molecule, where now system A is the He 
nucleus plus its electronic cloud on the left, and system B is the one on the right. We can take 
the two degenerate densities n2,i(r) and ni^2(r), 

n, i(r) = I ^K^^";) ^ ' ni ^(r) = | ^"^t^") " ^ 5?^ (12) 

'^^ I f^He+(r) re^B ' ^ ' 1 nHe(r) r e ^ ^ 

and form any ensemble of the two. We should always get the same energy. In particular, we can 
choose the symmetric one, ?^3/2,3/2(^) = |'^2,i(r) + ^'i-i,2(r), which yields | electrons on A and 
I electrons on B. Since the two degenerate wavefunctions corresponding to ?ii^2(r) and ra2,i(r) 
have zero overlap, we must have, for any < u; < 1, 

E[w 712,1 + {I -w) ni,2] = w E[n2,i] + {I - w) E[ni^2]- (13) 

If we now only look at the region pertaining to system A, we find |11|, [9] 



Ea[{1 - w) nHe+ +wnue] = (1 - w)-E^[nHe+] +wEA[n}ie], 



(14) 
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Figure 1. The non-interacting Kohn-Sham kinetic energy functional Tg, the exchange functional 
E'x; the correlation functional Ec, and their sum E^c for the ensemble n = {1 — w)ni + wn2 
formed by the A'^ = 1 and N = 2 densities of the Hooke's atom series, as a function of the weight 
w. All quantities are exact, calculated from the analytical solutions given by Taut [TT], and are 
reported in Hartree atomic units. 



which is the usual result for the energy of a fractional number of electrons [11^ [T2l I13j . 
corresponding to the well-known requirement for density functionals as a function of w, 

Ts[{l - w) UN + wnN-\-i\ + £^h[(1 - w) UN + wnN+i\ + £^xc[(l - w)nN if njv+i] = 

(1 - w) {Ts[nM] + EylIum] + E^cVin]) + w {Ts[nN+i\ + EuiuN+i] + E^c[nN+i\) , (15) 

where and 

?7.jv+i are the densities of the physical system (same Vne ) with iV and TV 1 
electrons, respectively. Equation psp is formally equivalent to Eq. (|10|) . We have only two 
densities n^r and htv+i, non-degenerate on the same region A, coming from the degeneracy 
arising in systems composed of many sub-systems [H]. 

Again, we may wonder whether Eq. p5|) can ever be attained by approximate functionals. 
To illustrate the complicated interplay between Tg and E^c as a function of w, we consider 
here the simple cases of the Hooke's atom series (two interacting electrons in an harmonic 
external potential, fne(r) = ^kr"^), and of the He isoelectronic series (two interacting electrons 
with t'ne(r) = The Hooke's atom series has a set of analytical solutions for specific 

values of the spring constant k [T^. Thus, we can calculate the exact densities and energies for 
A^ = 1 and N = 2 electrons, and, by inversion (see [H]), all the exact Kohn-Sham functionals 
corresponding to the ensemble density {1 — w)ni + w n2- In Figured] we report, as a function of 
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Figure 2. The non-interacting Kohn-Sham kinetic energy functional Tg, the exchange functional 
E'x; the correlation functional Ec, and their sum E^c for the ensemble n = {1 — w)ni + wn2 
formed by the N = 1 and N = 2 densities of the He atom series, as a function of the weight w. 
All quantities are very accurate, calculated from the variational wavefunctions of Ref. [19] (see 
also [20] and [2T]), and are reported in Hartree atomic units. 



w, the exact Ts[{l — w) ni + w 712] and E^c [(1 — z/;) ni + w n2] (with its two components, exchange 
and correlation, separately) for three different values of A:. As A: decreases, the system becomes 
more and more correlated. The non-interacting kinetic energy Tg as a function of w changes 
from being almost linear in the less correlated case k = j, to displaying a minimum in the more 

correlated ca.e t = 4 . 0.0003. The tunotronaj E„ mu=t compensate the qnadrafc 

behavior of En with w, 

Eu[{l-w)ni+wn2] = {1-wf Eu[ni] + w^ En[n2]+w{l-w) j dv j dv' ^^^^^^'^^^^ \ (16) 

as well as the non-linear behavior of Tg for correlated systems. 

A very similar trend of the functionals is observed for the He isoelectronic series, as shown 
in Fig. [21 In this case we have used an improved version [20] of the accurate variational 
wavefunctions and energies of Ref. [19] (see also [21]) to calculate the exact functionals. Notice 
in particular the behavior of Tg for the more correlated case, Z = \. 

A closed shell interacting system of two-electrons is weakly correlated when n2(r) ~ 2ni(r). 
In this case, the ensemble density with the corresponding = 1 system is simply given by 
wni + {l — w)n2 ~ {l + w)ni, and Tg[{l — w) ni + w n2] ~ {\ — w)Tg [rii] + wTg [n2] ■ This is the 



case, e.g., of the Ne*^ system of Fig. [21 for which Tg is almost hnear. However, the deviation of 
the Hartree term from hnearity, 

£'h[(1 - w)ni + wn2] - (1 - w)i?H["-i] - wEY{[n2] = 

w (l—w) cj r r / [n2(r)— ni(r)l[n2(r')— rai(r')l \ I 1 

= 4 — - J dr J dr' '-^^ — \r-r'\ — -^-^ = -w [1 - w) Eu_[n2 - ni\ (17) 

is significantly different from zero when n2 ~ 2ni, while it could become smaller when n2 is much 
more diffuse than 2 rii. For example, for an ensemble of H and H~ we have E}i[n2 — ni] = 0.1203 
Hartree, while Sni'i-i] = 0.3125 Hartree. That is, if the N = 2 system were less correlated so 
that 77-2 — ni were equal (or close) to ni, E'h would be further from the linear behavior. So (in 
the simple systems considered here) when Tg is closer to linearity E-n can be further from it, 
and viceversa. 

While the relative deviation of Tg from linearity decreases as the N = 2 system becomes less 
correlated (as shown in Figs. [T]and[2]), the maximum absolute value of Ts[{l — w)ni + wn2] — 
(1 — w)Ts[ni] — wTs[n2] is always of the same order of magnitude, i.e. ~ 2 mH for the Hooke's 
series and ~ 15 mHartree for the He series. 

5. Conclusions 

Size-consistency in DFT is often taken for granted with approximate functionals of the form of 
Eq. ([3]), because the density is believed to be an intensive quantity. However, the density is not 
intensive in the presence of degeneracy. An attempt to build a correct description of ensembles 
in DFT seems in order |22j . but using Eqs. (jlOp and ()15p is probably not the best starting point. 
The alternative path of building ensembles of Kohn-Sham systems (e.g., using different Kohn- 
Sham potentials for each symmetry of the degenerate system [ID] ) deserves further investigation. 
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